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Abstract 

Frozen orbits of the Hill problem are determined in the double averaged problem, where short 
and long period terms are removed by means of Lie transforms. The computation of initial con- 
ditions of corresponding quasi periodic solutions in the non-averaged problem is straightforward 
for the perturbation method used provides the explicit equations of the transformation that con- 
nects the averaged and non-averaged models. A fourth order analytical theory reveals necessary 
for the accurate computation of quasi periodic, frozen orbits. 


Introduction 

Besides its original application to the motion of the Moon [1], the Hill problem provides a good 
approximation to the real dynamics of a variety of systems, encompassing the motion of comets, 
natural and artificial satellites, distant moons of asteroids, or dynamical astronomy applications 
[2, 3, 4]. Specifically, the Hill model and its variations [5, 6, 7, 8, 9] are useful for describing motion 
about planetary satellites. Besides, the Hill problem is an invariant model that does not depend on 
any parameter, thus, giving broad generality to the results, whose application to different systems 
becomes a simple matter of scaling. 

A classical result shows that low eccentricity orbits around a primary are unstable for moderate 
and high inclinations due to third-body perturbations [10]. Almost circular orbits close to the 
central body remain with low eccentricity in the long term only when the mutual inclination with 
the perturbing body is less than the critical inclination of the third-body perturbations / = 39 . 2 ° 
(see [11] and references therein). Since low eccentricity, high inclination orbits are precisely the 
candidate orbits for science missions around natural satellites, a good understanding of the unstable 
dynamics of the Hill problem is required. 

The study of the long term dynamics is usually done in the double averaged problem. After 
removing the short and long period terms, the problem is of one degree of freedom in the eccentricity 
and the argument of the periapsis. As the double averaged problem is integrable, the phase space 
is made of closed curves and equilibria. The latter are orbits that, in average, have almost constant 
eccentricity and fixed argument of the periapsis, and are known as frozen orbits. 
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To each trajectory of the double averaged problem it corresponds a torus of quasi periodic 
solutions in the non-averaged problem. The accurate computation of initial conditions on the torus 
requires to recover the short and long period effects that were lost in the averaging. This is normally 
done by trial and error by making iterative corrections on the orbital elements, although other 
procedures can be applied [12]. 

Our analytical theory is computed with Deprit’s perturbation technique [13]. The procedure is 
systematic and has the advantage of providing the explicit transformation equations that connect 
the averaged analysis with proper initial conditions of the non-averaged problem. A second order 
truncation of the Hamiltonian shows that there are no degenerate equilibria and, therefore, is enough 
to give the qualitative description of the reduced system. However, the second order truncation 
introduces a symmetry of direct and retrograde orbits that is not part of the original problem, and 
a third order truncation is required to manifest the non-symmetries of the problem. 

While, in general, the third order theory provides good results in the computation of quasi 
periodic, frozen orbits, its solutions are slightly affected of long period oscillations. This fact may 
adversely affect the long term evolution of the frozen orbits and becomes apparent in the computation 
of science orbits about planetary satellites, a case in which small perturbations are enough for the 
unstable dynamics to defrost the argument of the periapsis. Then, the orbit immediately migrates 
along the unstable manifold with an exponential increase in the eccentricity. 

We find that a higher order truncation is desirable if one wants to use the analytical theory 
for computing accurate initial conditions of frozen orbits. The computation of the fourth order 
truncation removes almost all adverse effects from the quasi periodic solutions, and shows a high 
degree of agreement between the averaged and non-averaged models even in the case of unstable 
orbits. 

Whereas the third-body perturbation is the most important effect in destabilizing science orbits 
around planetary satellites, the impact of the non-sphericity of the central body may be taken into 
account. Previous research including both effects has limited up to third order theories (see [14] and 
references therein), but from the conclusions of this paper it may worth to develop a higher order 
theory including the inhomogeneities of the satellite’s gravitational potential. 


Double averaged Hill problem to the fourth order 

The equations of motion of the Hill problem are derived from the Hamiltonian 

H={\/2){X-X)-u-{xy.X) + W{x), W = (w 2 /2)(r 2 - 3 .t 2 ) - g/r (1) 

where, in the standard coordinate system of Hill’s model, x = (. x , y, z ) is the position vector, X = 
( X,Y,Z ) is the vector of conjugate momenta, r = ||x||, and both the rotation rate of the system 
to = 1 1 a; 1 1 and the gravitational parameter /i of the primary are set to 1 in appropriate units. 

The problem is of three degrees of freedom, yet admitting the Jacobi constant H = —C /2. 
Despite its non-integrability, approximate solutions that explain the long term dynamics can be 
found by perturbation methods. Close to the central body the Hill problem writes as the perturbed 
two body problem 

H = (1/2) (X 2 + Y 2 + Z 2 ) - (1/r) - e(xY - y X) + (e 2 /2) (r 2 - 3x 2 ), (2) 

where e is a formal parameter used to manifest the order of each perturbation of the Keplerian 
motion. Thus, the Coriolis term is a first order effect and the third body perturbation appears at 
the second order. 

To apply perturbation theory, we formulate the problem in Delaunay variables {£, g , h, L , G, H), 
where £ is the mean anomaly, g is the argument of the periapsis, h the argument of the node in 
the rotating frame, L = y/JTa is the Delaunay action, G = L \/l — e 2 is the modulus of the angular 
momentum vector, H = Geos / is its polar component, and a, e, /, are usual orbital elements: 
semimajor axis, eccentricity and inclination. 
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Our theory is based on the use of Lie transforms as described by Deprit [13, 15]. It has the advan- 
tage of connecting the averaged and original problems through explicit transformations equations. 
After removing the short and long period terms we get the transformed Hamiltonian 

1C = K 0 ,o + £ K 0 ,i + (£ 2 /2) Ko, 2 + (e 3 /6) K 0>3 + (e 4 /24) K 0 , 4 (3) 


where e = L 3 , 

K 0 ,o = -1/(2 A 2 ), 

Ko,i = Ao,o 2<r 

/v 0 , 2 = A'o.o (1/4) [(2 + 3e 2 ) (2 - 3s 2 ) + 15e 2 s 2 cos 2 g] , 

K 0 , 3 = ATo, 0 (27/32) a [2s 2 + (50 - 17s 2 ) e 2 + 15e 2 s 2 cos 2 g\ , 

K oa = K 0 ,o (-3/512) {3285s 4 e 4 cos 4 g - 12s 2 [3996 - 2940s 2 - (4582 - 4035s 2 ) e 2 ] e 2 cos 2 g, 

+8(784 - 708s 2 - 9s 4 ) - 144(926 - 941s 2 + 244s 4 ) e 2 + 9(10728 - 15208s 2 + 5007s 4 ) e 4 } , 


and cr = H / L = ci], r) = -\/l — e 2 is the eccentricity function, and we use the abbreviations s = sin/, 
c = cos /. Details on the perturbation method and expressions to compute the transformation 
equations of the averaging are given in the Appendix. 

The double averaged Hamiltonian, Eq. (3), depends neither on the mean anomaly nor on the 
argument of the node. Therefore, the corresponding conjugate momenta, L and H , are integrals of 
the reduced problem and Hamiltonian (3) represents a one degree of freedom problem in g and G. The 
equations of motion are computed from the Hamilton equations dg/dt = dIC/dG , dG/dt = —dIC/dg, 


dg 

dt 


dG 

dt 


[5c 2 — rf — 5 (c 2 — rj 2 ) cos 2 g] + - [5c 2 + llrj 2 — 5 (c 2 — rj 2 ) cos 2 g] 

(jr 4Gr 

3e 2 r 

+ — — { 2113c 2 - 3285c 4 + (3915 + 9165c 4 ) rj 2 + (1581 + 7791c 2 ) r ] 4 
128G L 

-4 [802c 2 - 1095c 4 + (19 + 2565c 4 ) if - (547 + 1744c 2 ) ry 4 ] cos2 ff 
+1095e 2 s 2 (c 2 -ry 2 )cos4.g|, 


— ^e 2 s 2 |5 (8 + 9ecr) sin 2 g 
e 2 / 

+ — {2 [509 - 1095c 2 + (547 + 4035c 2 ) rf] sin 2g 


1095 e 2 s 2 sin 4 g'j |. 


( 4 ) 


( 5 ) 


Once integrated g and G for given initial conditions, the secular variations of £ and h are computed 
from simple quadratures derived from the Hamilton equations dh/dt = dK./dH , d£/dt = dlC/dL 


h — ho -\- 


d_ 

dH 


IC(g(t),G(t);H,L)dt, 


£ = 


d_ 

dL 


lC(g(t),G(ty,H, L) dt. 


Qualitative dynamics 

The flow can be integrated from the differential equations above, Eqs. (4)-(5). However, since the 
system defined through Eqs. (4) and (5) is integrable the flow is made of closed curves and equilibria, 
and it can be represented by contour plots of Hamiltonian (3). Thus, for given values of the dynamic 
parameters L and H — or £ and a - we can plot the flow in different maps that are function of g , G. 
Figure 1 shows an example in semi-equinoctial elements (e cos g,e sing), where we note a hyperbolic 
point corresponding to an unstable circular orbit, and two elliptic points corresponding to two stable 
elliptic orbits with e = 0.3 and periapsis at g = ±7t/2 respectively. 
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L=0 . 295105 , H=0. 215826 



-1 -0.5 0 0.5 1 

e*cos (g) 


Figure 1: Flow in the doubly reduced phase space. 


Delaunay variables are singular for zero eccentricity orbits, where the argument of the periapsis 
is not defined. Then, it is common to study the reduced phase space in the variables introduced by 
Coffey, Deprit, and Deprit, [16, 8] 

Xi = r/escosg, % 2 = V e s sin g, X 3 = V 2 - (1/2) (1 + cr 2 ), 

that define the surface of a sphere of radius R = (1/2) (1 — cr 2 ). 1 Then, after dropping constant 
terms and scaling, Hamiltonian (3) writes 


K. = -127? 2 - + ^ecr ( 25- 24?y 2 -cr 2 - 15 


V 


e 2 - 
64 L 


3815 + 9528cr 2 + 9cr 4 


( 6 ) 


-6 (343 + 1709cr 2 ) rf - 1824r/ 4 + 6 (293 - 821^ - 1470cr 2 ) - 3285 


xi 


V 


V 


The flow on the sphere is obtained from Liouville equations Xi = {Xij K - } > * = 1,2, 3, where the dot 
means derivative in the new time scale. Then, 


Xi = 


X2 = 


X3 = 


I677 


X 2 <64 (5 - 8 rf + 5ct 2 ) + 72e cr (5 - 2rf + 5cr 2 ) + — 3815 - 1824? ? 4 + 9528cr 2 (7) 


64 L 


+9ct 4 - Qrf (343 + 1709ct 2 ) + 6 (293 - 82 I 77 2 - 1470cr 2 ) (x 2 /*?) 2 - 3285 (X 2 /*?) 4 ] } , 

— — Xi ( 608e 2 rf A + rf \ 128 + 576e a + e 2 (343 + 1709ct 2 )1 
lOT] l L J 

- [320 + 360e cr - e 2 (293 - 1470ct 2 )] (x 2 /*?) 2 - 1095 e 2 (x 2 /??) 4 }, 

^ Xi X 2 {320 + 360e cr - e 2 [293 - 1470cr 2 - 82 I 77 2 - 1095 (x' 2 A?) 2 ] } , 


( 8 ) 


(9) 


with the constrain xi Xi + X 2 X 2 + X 3 X 3 = 0. 

x The sphere representation misses the case G = H = 0, irrelevant in astrodynamics. Full details including 
rectilinear solution will be given in a future paper. 
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Equations (7)-(9) show that circular orbits (xi = X'i = 0, Xa = R, the “north” pole of the 
sphere) are always equilibria. Equations (8) and (9) vanish when \i = 0, / 0, but not Eq. (7), 

which vanishes only when 

1095e 2 cr 4 - a 2 [320 + 36(k a + e 2 (802 + 2565a 2 )] if (10) 

+ [l92-216£a-e 2 (362 - 35a 2 )] rf - &le 2 rf = 0. 

Equation (10) is a polynomial equation of degree 8 in r /, therefore admitting eight roots. Note that, 
for the accepted values of e <C 1, Eq. (10) is of the form Af — A\ x + A\ x 3 — A\ x 4 = 0 that admits 
a maximum of two real roots, according to Descartes’ rule. 

The real roots of Eq. (10) that verify the dynamical constrain |a| < r/ < 1 are also equilibria. 
The root i] = 1 marks a change in the number of equilibria or “bifurcation” (rj > 1 could be a root 
but not an equilibrium). Then, the number of equilibria changes along the line 

„ -27a - 45a 3 ± V5076 + 1473a 2 + 4730a 4 - 27375a 6 

e = 4 s (11) 

423 + 767a 2 + 1470a 4 v ' 

obtained making 77 = 1 in Eq. (10), that establishes a relation between the dynamical parameters 
e and a corresponding to bifurcations of circular orbits. Figure 2 shows that this line defines two 
regions in the parameters plane with different number of equilibria in phase space. Circular orbits 


Bifurcation line of circular orbits 



Figure 2: Regions in the parameters plane with different numbers of equilibria. 

in the outside region of the curve are stable. When crossing the line given by Eq. (11) the number 
of real roots of Eq. (10) with dynamical sense increases in a pitchfork bifurcation: circular orbits 
change to instability and two stable elliptic orbits appear with periapsis respectively at g = ±7r/2, 
as in the example of Figure 1. 

Note that the curve given by Eq. (11) notably modifies the classical inclination limit cos 2 I > 3/5 
for circular orbits’ stability. However, we can not extend the practical application of the analytical 
theory to any value of e. It is common to limit the validity of the Hill problem approximation to 
one third of the Hill radius rn = 3 -1 / 3 . Then e < (rn/3) 3 ^ 2 = 1/9, including most of the planetary 
satellites of interest. Figure 3 shows the bifurcation lines of circular orbits in the validity region of the 
parameters plane with the e values corresponding to low altitude orbits around different planetary 
satellites highlighted. 


Frozen orbits computation 

Hill’s case of close orbits to the smaller primary is a simplification of the restricted three-body 
problem, which in turn is a simplification of real models. Therefore, the final goal of our theory is 
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Figure 3: Bifurcation lines in the parameter plane. 


not the generation of ephemerides but to help in mission designing for artificial satellite missions 
about planetary satellites, where frozen orbits are of major interest. 

For given values of the parameters e, a , determined by the mission, a number of frozen orbits may 
exist. A circular frozen orbit, either stable or unstable, exists always and the computation of real 
roots \a\ < g < 1 of Eq. (10), if any, will provide the eccentricities of the stable elliptic solutions with 
frozen periapsis at g = ±7r/2. To each equilibrium of the doubly reduced phase space corresponds 
a torus of quasi periodic solutions in the original, non-averaged model. Below we present several 
examples that justify the effort in computing a fourth order theory to reach the quasi periodicity 
condition in the Hill problem. 

Elliptic frozen orbits 

We choose e = 0.0470573, a = 0.422618. If we first try the classical double averaged solution, 
the Hamiltonian Eq. (3) is simplified to ATo,o + £ Ko,i ± (er 2 /2) Kq, 2 , and the existence of elliptic 
frozen orbits reduces to a 2 < 3/5, g = ±7r/2. The eccentricity of the elliptic frozen solutions is 
then computed from = (5cr 2 /3) 1 ^ 4 — obtained by neglecting terms in e in Eq. (10). Thus, for the 
given values of e and <r, and taking into account that we are free to choose the initial values of the 
averaged angles f, h , we get the orbital elements of the first row of Table 1. The left column of Fig. 
4 shows the long term evolution of the instantaneous orbital elements for this case, in which we find 
long period oscillations of more than four degrees in inclination, more than fifteen in the argument 
of periapsis, and a variation of ±0.06 in the eccentricity. 

When computing a second order theory with the Lie-Deprit perturbation method we arrive 
exactly at the classical Hamiltonian and bifurcation condition. However, we have available the 
transformation equations to recover the short and long period effects, although up to the first order 
only. After undoing the transformation equations we find the orbital elements of the second row 
of Table 1 where we see that all the elements remain unchanged except for the eccentricity and 
inclination. The long term evolution of these elements is presented in the right column of Fig. 
4, in which we notice a significant reduction in the amplitude of long period oscillations: 2.5° in 
inclination, around 10° in the argument of the periapsis, and ±0.04 in eccentricity. 

The results of the third and fourth order theories are presented in the last two rows of Table 1 and 
in Fig. 5. The higher order corrections drive slight enlargements in the semimajor axis. While both 
higher order theories produce impressive improvements, we note a residual long period oscillation 
in the elements computed from the third order theory (left column of Fig. 5). On the contrary, the 
orbital elements of the frozen orbit computed with the fourth order theory are almost free from long 
period oscillations and mainly show the short period oscillations typical of quasi periodic orbits. 
Furthermore, after applying standard differential correction to the initial conditions computed from 
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the fourth order theory we easily converge to an orbit of the non-averaged problem, presented in 
Fig. 8, that is periodic after 99 cycles. 

Table 1: Initial orbital elements of an elliptic frozen orbit for e = 0.0470573, cr = 0.422618. 


Theory 

a (Hill units) 

e 

/ (deg) 

9 (deg) 

h 

i 

Classical 

0.130342 

0.674094 

55.0995 

-90 

0 

0 

2nd order 

0.130342 

0.648065 

55.6915 

-90 

0 

0 

3rd order 

0.130515 

0.637316 

56.1798 

-90 

0 

0 

4th order 

0.130538 

0.634803 

56.2813 

-90 

0 

0 


Circular frozen orbits 

If we choose the same value for e but now a = 0.777146, frozen elliptic orbits do not exist any longer 
and the circular frozen orbit is stable. Both the third and fourth order theories provide good results, 
but, again, the third order theory results in small long period oscillations in the eccentricity whereas 
the fourth order theory leads to a quasi periodic orbit (see Fig. 6). 

For e = 0.0339919 and cr = 0.34202 the circular frozen orbit is unstable. Due to the instability, 
a long term propagation of the initial conditions from either the third or the fourth theory shows 
that the orbit escapes by the unstable manifold with exponential increase in the eccentricity. But, 
as Fig. 7 shows, the orbit remains frozen much more time when using the fourth order theory. A 
variety of tests performed on science orbits close to Galilean moons Europa and Callisto showed 
that the fourth order theory generally improves by 50% the lifetimes reached when using the third 
order theory. 

As in the previous case of a elliptic orbit, standard differential corrections converge to the periodic 
orbits presented in Fig. 9. 


Conclusions 

Frozen orbits computation is a useful procedure in mission designing for artificial satellites. After 
locating the frozen orbit of interest in a double averaged problem, usual procedures for computing 
initial conditions of frozen orbits resort to trial and error interactive corrections, or require involved 
computations. However, the explicit transformation equations between averaged and non-averaged 
models can be obtained with analytical theories based on the Lie-Deprit perturbation method, which 
makes the frozen orbits computations straightforward. 

Accurate computations of the initial conditions of frozen, quasi periodic orbits can be reached 
with higher order analytical theories. This way of proceeding should not be undervalued in the 
computation of science orbits around planetary satellites, a case in which third-body perturbations 
induce unstable dynamics. 

Higher order analytical theories are a common tool for computing ephemeris among the celestial 
mechanics community. They are usually developed with specific purpose, sophisticated algebraic 
manipulators. However, the impressive performances of modern computers and software allow us 
to build our analytical theory with commercial, general-purpose manipulators, a fact that may 
challenge aerospace engineers to use the safe, well known techniques advocated in this paper. 
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Appendix 

Let T : ( x,X ) — > ( x',X '), where x are coordinates and X their conjugate momenta, be a Lie 
transform from “new” (primes) to “old” variables. If W = ^^(e’/i!) W»+ \{x,X) is its generating 



function expanded as a power series in a small parameter e, a function T = JT(e*/i!) Fj,o(cc, -X") 
can be expressed in the new variables as the power series (T : T) = /i\) Fo t i(x' , X 1 ) whose 

coefficients are computed from the recurrence 

Fid=Fi+U-i+ E (k) {Fk,j-i;W i+1 - k } (12) 

0 <k<i ' ' 


where {Fkj-i\Wi + i-k} = V x Fk,j - 1 • Vy W i+1 - k - Vjftj-i • Wi + i_ fc , is the Poisson bracket. 
Conversely, the coefficients Wj+i of the generating function can be computed step by step from Eq. 

(12) once corresponding terms F o i of the transformed function are chosen as desired. In perturbation 
theory it is common to chose the #o,i as an averaged expression over some variable, but it is not the 
unique possibility [17]. Full details can be found in the literature [18, 19]. 

To average the short period effects we write Hamiltonian (2) in Delaunay variables as 

n = # 0 ,o + e Hi t o + (e 2 /2) # 2 , o + (e 3 /6) # 3,0 + (e 4 /24) #4, o (13) 

where #o,o = — 1/ (2i 2 ) , #i , 0 = — #, FI 2j0 = r 2 {l — 3 [cos(/ + g) cos h — csin(/ + g) sin h] 2 }, and 
# 3 ,o = F [ 4 0 = 0. Note that the true anomaly / is an implicit function of t. 

Since the radius r never appears in denominators, it results convenient to express Hamiltonian 

(13) as a function of the elliptic — instead of the true — anomaly u by using the ellipse relations 
r sin / = 77 a sin u, r cos / = a (cos u — e), r = a (1 — ecos it). 

After applying the Delaunay normalization [20] up to the fourth order in the Hamiltonian, we 
get 

H' = H 0 ,o + e Hop + (e 2 /2) H 0 , 2 + (e 3 / 6 ) H 0i 3 + (e 4 /24) H 0A (14) 

where, omitting primes, 

-1/(2 L 2 ) 

H 0 fi e2cq 

tf 0 , 0 (e 2 / 8 ) { ( 4 + 6e 2 ) (2 - 3s 2 + 3s 2 cos 2 h) 

+15 e 2 [2s 2 cos 2g + (1 — c) 2 cos(2g — 2 h) + (1 + c) 2 cos(2 g + 2h)] } 

# 0,0 (45e 3 /8) e 2 77 [(1 + c ) 2 cos(2 g + 2 h) — (1 — c ) 2 cos(2 g — 2 ft.)] 

# 0,0 (— 3e 4 /512) 1 16 (47 + 282c 2 + 63c 4 ) - 144 (227 + 90c 2 + 59c 4 ) e 2 
— 18 (227 + 610c 2 - 701c 4 ) e 4 - 24s 2 [558 + 270c 2 + (109 - 555c 2 ) e 2 ] e 2 cos 2 g 

+24s 2 [216 + 56c 2 - 8 (161 + 59c 2 ) e 2 - (11 - 701c 2 ) e 4 ] cos 2 h 
-48 (1 + c ) 2 [338 - 90c + 90c 2 - (91 - 185c + 185c 2 ) e 2 ] e 2 cos(2 g + 2 h) 

-48 (1 - c ) 2 [338 + 90c + 90c 2 - (91 + 185c + 185c 2 ) e 2 ] e 2 cos(2 ff - 2 h) 

+ 6 s 4 (56 - 472e 2 + 701e 4 ) cos 4 h + 1710s 4 e 4 cos 4 g 
—60s 2 (18 — 37e 2 ) e 2 [(1 + c ) 2 cos(2 g + Ah) + (1 — c ) 2 cos(2 g — Ah)] 

+1140s 2 e 4 [(1 + c ) 2 cos(4g + 2 h) + (1 — c ) 2 cos(4g — 2ft,)] 

+285 e 4 [(1 + c ) 4 cos(4g + 4ft) + (1 — c ) 4 cos(4t; — 4ft)] j 

The generating function of the transformation is W = W 2 + (1/2) W 3 , where 

W 2 = L(£ 2 /192){4(2-3s 2 ) [3e(5 + 377 2 )^ 1 ,o,o- 9 e 2 5 2 ,o,o + e 3 5 3 ,o,o] 

+ 6 s“ e [3 (5 + 3 t; 2 ) ( 51 , 0,2 + 5i,o- 2 ) — 9e (^ 2 , 0,2 + 5 2 ,o- 2 ) + e 2 (^ 3 , 0,2 + ^ 3 , 0 - 2 )] 

+ 6 s 2 (1 + g ) 2 [15e 5 1>2 , 0 - (9 - 677 ) 5 2 , 2 , 0 + e 5 3 , 2 , 0 ] 

+ 6 s 2 (1 - 77) 2 [15e 5i _ 2)0 - (9 + 677 ) 5 2 ,_ 2 , 0 + e 5 3 ,_ 2 , 0 ] 


ffo,o = 

#0.1 = 

Ho, 2 = 

Ho, 3 = 
#0,4 = 
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+3 (1 + c) 2 (1 + rj ) 2 [15e 51,2,2 — (9 — 6 if) 82 , 2,2 + e 53 , 2 , 2 ] 

+3 (1 - c) 2 (1 + rj ) 2 [15e 5i, 2,-2 - (9 - 6??) S 2 , 2,-2 + e 5 3 , 2 - 2 ] 

+3 (1 — c ) 2 (1 — rj ) 2 [15e 5 i,_ 2,2 — (9 + 677 ) S 2 ,~ 2 , 2 + e 5 3 ._ 2 , 2 ] 

+3 (1 + c ) 2 (1 — rj)~ [15e 5 i,_2,-2 — (9 + 6 ij) S 2 - 2,- -2 + e <S' 3 ,_ 2 ,— 2 ] } 

W 3 = L (e 3 /256) {72e s 2 (13 + 3r; 2 ) [5 i, 0 ,2 - Si, 0 ,- 2 ] - 24e 2 s 2 (17 + 4, ? 2 ) [5 2 ,o,2 - 5 2 ,o,- 2 ] 
+88e 3 s 2 [5 3 ,o,2 — 5 3 ,o,— 2 ] — 6e 4 s 2 [ 54 , 0,2 — 54, 0 ,- 2 ] 

+36 e (1 + 77 ) (13 + ry + 8 ? 7 2 ) [(1 + c ) 2 5 i, 2,2 ~ (1 — c ) 2 5 i, 2 ,— 2 ] 

+36 e (1 - rj) (13 - 77 + 8 rj 2 ) [(1 - c ) 2 5i ,_ 2 ,2 - (1 + c ) 2 Si ,- 2 ,- 2 ] 

-12 (1 + 77) 2 (17 - 677 - 877 2 ) [(1 + c ) 2 5 2 , 2 ,2 - (1 - c ) 2 5 2 , 2 ,- 2 ] 

-12 (1 - T ?) 2 (17 + 677 - 8 ? 7 2 ) [(1 - c ) 2 5 2 , _ 2,2 - (1 + c ) 2 5 2 ,- 2 ,- 2 ] 

+4 (1 + 77) 2 e (11 - 677 ) [(1 + c ) 2 5 3 , 2 ,2 - (1 - c ) 2 5 3 , 2 ,- 2 ] 

+4 (1 - 77) 2 e (11 + 677 ) [(1 - c ) 2 5 3 ,_ 2,2 - (1 + c ) 2 5 3 ,_ 2 ,- 2 ] 

—3 (1 + rj ) 2 e 2 [(1 + c) 2 54 , 2,2 — (1 — c) 2 54, 2 , - 2 ] 

-3 (1 - 77) 2 e 2 [(1 - c ) 2 5 4 ,_ 2,2 - (1 + c ) 2 5 4 ,_ 2 ,— 2 ] } 

We shorten notation calling Sij,k = sin(i u + j g + kh). 

The Lie transform of generating function W can be applied to any function of Delaunay variables 
T = JT(e® /*!) Fi{l ' , g' ,h ' , L' ,G ' , H'). Since W\ = 0, up to the third order in the small parameter 
recurrence ( 12 ) gives 

T = F 0 + (e 2 /2) {F 0 ; W 2 } + (e 3 /6) {F 0 ; W 3 } 

Specifically, this applies to the transformation equations of the Delaunay variables themselves, where 
+0 £ (P , g'i h’ , L', G' , H') and 5) = 0 for i > 0. 

A new application of the recurrence Eq. ( 12 ) to the Hamiltonian K, = 5 Zo<i <4 ( e V*0 0 , where 

Ki,o = Ho,i of Eq. (14), allows to eliminate the node up to the fourth order, obtaining the double 
averaged Hamiltonian Eq. (3). Note that A' 0,4 corrects previous results in [21]. The generating 
function of the transformation isV’ = V - i + eV 2 + (e 2 / 2) V 3 , where, omitting double primes, 

V\ = L (3/64) [(4 + 6 e 2 ) s 2 sin 2h + 5 (1 + c ) 2 e 2 sin(2g + 2 h) — 5 (1 — c ) 2 e 2 sin( 2(7 — 2 h)] 

V 2 = L (-3/128) 77 [ 6 c (2 - 17e 2 ) s 2 sin 2 h + 5 (2 - 9c) (1 + c ) 2 e 2 sin(2 5 + 2 h) 

+5 (1 — c ) 2 (2 + 9c) e 2 sin(2 g — 2 h) ] 

V 3 = L (-9/32768) { 16s 2 [456 - 104c 2 - 8 (193 + 754c 2 ) e 2 + (47 + 7831c 2 ) e 4 ] sin 2 h 
+2 s 4 (232 + 416e 2 - 1803e 4 ) sin Ah 

-32 (1 + c ) 2 e 2 [2 (323 - 285c + 780c 2 ) - (527 - 1135c + 2125c 2 ) e 2 ] sin(2 ff + 2 h) 

+32 (1 - c ) 2 e 2 [2 (323 + 285c + 780c 2 ) - (527 + 1135c + 2125c 2 ) e 2 ] sin(2 ff - 2h) 
+220s 2 e 2 (4 — lie 2 ) [(1 + c ) 2 sin( 2<7 + 4 h) — (1 — c ) 2 sin( 2<7 — 4 h)] 

+4520s 2 e 4 [(1 + c ) 2 sin(4g + 2 h) — (1 — c ) 2 sin(4 g — 2 h)] 

— 385e 4 [(1 + c ) 4 sin(4g + 4 h) — (1 — c ) 4 sin(4g — 4/i)] } 

The new Lie transform of generating function V can be applied to any function of Delaunay vari- 
ables, and, specifically, to the Delaunay variables themselves. For any £" € [t" ,g" 7 h" ,L” ,G" ,H”) 
the transformation equations of the Lie transform are computed, up to the third order, from 

+ £ <5i + (e 2 /2) 82 + (e 3 /6) <5 3 , 

where 

Si = U"; Vi}, 82 = v 2 j + {^; Vi}, 5 3 = {£"; ^3} + {{?"; ^2}; 4 i} + {<y i; v 2 } + Vi} 
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Figure 8: Elliptic, stable, periodic orbit after 99 cycles. 
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Figure 9: Left: Near circular, stable, periodic orbit after 63 cycles. Right: Near circular, unstable, 
periodic orbit after 60 cycles. 
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